This invention relates to the field of image processing, and in particular to the calculation of local variance in coherent images such as ultrasound images used in medical and non-destructive testing applications.
Ultrasound imaging systems for medical diagnostics scan a subject with ultrasound beams in a sector or rectangular shape and provide a two-dimensional display of the scanned field. The imaging system includes a transducer that emits and focuses the ultrasound beams along certain linear directions (scan lines). The emitted beams are refracted, absorbed, dispersed and reflected within the target subject. Those beams that are reflected back to the transducer are converted to analog electrical signals. These signals undergo a number of processing steps such as filtering, amplification, analog to digital conversion, beamforming, and demodulation. These processing steps produce a two-dimensional array of digital data representative of the acoustic properties of the material at various points in the scanning field. This array is then processed by image processing circuitry. The image processing circuitry removes imperfections and improves the desired features of the image. The image processing circuitry also encodes the beam-formed formed data into picture elements (pixels). The pixels are then converted to an analog signal for output to a viewing device, such as a CRT monitor. The resulting B-mode image (brightness image) is a gray scale image in which the range of darker-to-lighter gray-scale shades correspond to increasing strength of the reflected signal (echo).
An ultrasound image contains innumerable imperfections including electronic and acoustic noise, misplaced information due to refraction and scatter of ultrasound beams and, under some circumstances, gaps in the image which contain no genuine ultrasound information. Several types of image processing methods have been devised to manipulate the digital data, permitting improvement in the final displayed image. In some image processing methodsxe2x80x94such as filtering for speckle reduction in ultrasound imagesxe2x80x94it is necessary to calculate various local statistical characteristics at each datum in the image. Local statistical characteristics are those statistical characteristics calculated using the data contained in a rectangular window around the sample of interest, as opposed to global characteristics which are calculated using all the data in the array.
Two of the most common local characteristics used in image processing applications are the mean and the variance, defined mathematically by the formulas:                               MEAN                      n            ,            m                          =                              1                          H              *              V                                ⁢                                    ∑                              p                =                0                                            V                -                1                                      ⁢                          xe2x80x83                        ⁢                                          ∑                                  q                  =                  0                                                  H                  -                  1                                            ⁢                              xe2x80x83                            ⁢                              S                                                      n                    +                    p                                    ,                                      m                    +                    q                                                                                                          (        1        )                                          VAR                      n            ,            m                          =                              1                          H              *              V                                ⁢                                    ∑                              p                =                0                                            V                -                1                                      ⁢                          xe2x80x83                        ⁢                                          ∑                                  q                  =                  0                                                  H                  -                  1                                            ⁢                              xe2x80x83                            ⁢                                                (                                                            S                                                                        n                          +                          p                                                ,                                                  m                          +                          q                                                                                      -                                          MEAN                                              n                        ,                        m                                                                              )                                2                                                                        (        2        )            
Where, as described for FIG. 1, xe2x80x9cHxe2x80x9d and xe2x80x9cVxe2x80x9d are the horizontal and vertical size (in samples) of the rectangular window in which the statistics are calculated, xe2x80x9cnxe2x80x9d and xe2x80x9cmxe2x80x9d are the line (vertical) and column (horizontal) indexes of the origin of the window relative to the origin of the array, xe2x80x9cpxe2x80x9d and xe2x80x9cqxe2x80x9d are the vertical and horizontal position indexes relative to the window origin, and Sij is the sample at position (ij).
Typically these characteristics are calculated sequentially for every sample Sij in the image, starting at the origin of the image (j=0,i=0) then moving the window one sample position at a time in the line (xe2x80x9cjxe2x80x9d ) direction until the local characteristics are calculated for every sample in the line. The window is then returned to the start of the lines (j=0) and moved down one position in the column direction (xe2x80x9cixe2x80x9d) to begin calculating the local characteristics for every sample in the next line. This pattern is repeated until the local characteristics are calculated for every sample in the image.
The characteristics have to be calculated in xe2x80x98real-timexe2x80x99. In other words, the characteristics have to be calculated at the same rate that the samples are generated in order to display typically several tens of images per second. Since the number of samples in each image is generally large (tens or hundreds of thousands), the real-time sample rate tends to be high, on the order of tens of millions per second. This leads to a large computational load for the apparatus used to compute the local characteristics. In order to keep the cost of the apparatus reasonably low, it is necessary to devise methods to calculate the characteristics with as few and as simple operations as possible.
Unfortunately, with the aforementioned algorithm, the direct calculation of the mean using formula (1) requires H*Vxe2x88x921 additions, and the direct calculation of the variance using formula (2) requires H*Vxe2x88x921 additions, H*V subtractions and H*V squaring operations; for example, if H=V=10, direct calculation of the mean would require 99 additions and the direct calculation of the variance would require 199 addition/subtractions and 100 squarings. Because these calculations have to be performed for each sample in the image, the direct calculation of mean and variance cannot be performed in real-time at a reasonable cost with the current technology.
Recursive sum algorithms are known and have been used for a long time for reducing the number of computations necessary to calculate the local mean. In these algorithms the window is moved over the array one sample position at a time and the sum of samples for the new window position is calculated from the sum of the previous position by subtracting the samples which drop out of the window and adding the new samples which enter the window. No calculations are needed for the samples belonging to both the old and the new window position. Therefore, only two additions and two subtractions per window position are performed for a 2-dimensional window of any sizexe2x80x94a considerable saving when compared with the 99 additions for the 10 by 10 window in example above.
However, it is not possible to directly implement the recursive sum algorithms to the calculation of local variance using formula (2). This is because, in formula (2), all terms in the summation   (            ∑              p        =        0                    V        -        1              ⁢          xe2x80x83        ⁢                  ∑                  q          =          0                          H          -          1                    ⁢              xe2x80x83            ⁢                        (                                    S                                                n                  +                  p                                ,                                  m                  +                  q                                                      -                          MEAN                              n                ,                m                                              )                2              )
change with every change in the window position. All the terms depend not only on the values of the samples in the window, but also on the mean, which changes with the window position.
The above-discussed and other drawbacks and deficiencies of the prior art are overcome or alleviated by the apparatus and method for the real-time calculation of local variance in images of the present invention. In an exemplary embodiment, a method of calculating local variance in an array of digital signals includes defining a window for moving within the array of digital signals, with the window having a size to include at least two digital signals. First, the window is moved to a first position where the window includes first and second digital signals, with the first and second digital signals indicative of first and second values. At this first position, the first and second values are added to obtain a first sum. The square of the first and second values are also added to obtain a second sum. A variance value is then determined for the first window position using the window size, the first sum, and the second sum. Next, the window is moved to a second position where the window includes the second digital signal and a third digital signal, with the second digital signal indicative of the second value and the third digital signal indicative of a third value. At the second position, the first sum and the third value are added, then the first value is subtracted to obtain a third sum. The second sum and the square of the third sum are added, then the square of the first sum is subtracted to obtain a fourth sum. Finally, a variance value for the second window is determined using the window size, the third sum and the fourth sum.
In an alternate embodiment, the method of calculating local variance in an array of digital signals includes defining a window for moving within an array of digital signals, with the window having a size to include at least two digital signals. A plurality of values indicated by a plurality of signals within the window are recursively summed to obtain a first window sum. The plurality of values within the window are squared, and the squares of the plurality of values are recursively summed to obtain a second window sum. The first and second window sums are then used to determine a variance value for the window position. The variance value is determined by, first, dividing the square of the first window sum by the size of the window to obtain a first interim value. Next, the first interim value is subtracted from the second window sum to obtain a second interim value. Finally, the second interim value is divided by the size of the window to obtain the variance value for the window position.